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1. Introduction 

Population dynamics is one of the natural areas of application of the theory of stochastic 
processes: Parameters of deterministic models are perturbed by random fluctuations in 
order to model the influence of an ever-changing environment consisting of a multitude 
of individual "agents". Analytical solutions are known for only a handful of such 
systems. A majority of physically important systems is accessible only via numerical 
simulations. It is, therefore, particularly important to thoroughly understand those few 
whose analytical properties are known and test how well the existing numerical methods 
can approximate them. When a rigorous mathematical statement can be made about a 
system described by a stochastic differential equation (SDE), one expects that it will be 
corroborated by numerical simulations or by an actual physical experiment, if the latter 
is possible: after all, an SDE is a mathematical idealization of a physical reality that is 
discrete by its very nature. If the results of the two approaches, continuous time SDE 
vs. discrete time numerical simulations, do not agree, there are reasons for a serious 
concern. 

One of the most surprising results, reported recently by Mao and co-workers in 
Ref. P, is the fact that a multiplicative noise can stabilize a system whose deterministic 
counterpart is divergent. The class of systems discussed in that Reference may comprise 
several interacting species but its simplest form, to which we will restrict ourselves within 
this paper, is described by a scalar Langevin equation 



where r > and is a Gaussian white noise (GWN), interpreted in the sense of 
Ito, with = and (^(^1)^(^2)) = 8(ti — t 2 ). If there is no noise, a = 0, the 

"— " sign correspond to the well-known logistic equation, while the "+" sign describes a 
"population" that diverges in a finite time. This latter case defies one of the ecological 
laws that a population should be self-limited 0, although Eq. ([I]) perhaps can be used 
to model explosive phenomena observed in some chemical, nuclear or stellar systems 
where a divergence is interpreted as a destruction of the original system. However, our 
primary purpose of discussing the "+" sign is to present the suprising effects that noise 
can have on a deterministically divergent system. 

It is rigorously shown in Ref. pQ| that for any nonzero noise, a 7^ 0, if the system Q 
starts from a positive initial value, it will for all times remain positive and bounded 
almost surely, regardless of the choice of the sign inside the brackets. This means that 
even a tiny addition of the noise can stabilize the system. We note that a proof of this 
statement runs in two stages: First it is proven that the system remains positive, and 
later this fact is used to prove that the system remains bounded. 

Apart from the formal proof, Ref. pQ provides several numerical examples. These 
examples, however, are much less convincing than the fine mathematical points they are 
supposed to illustrate. First, it is not clear whether the numerical trajectories presented 
are typical or handpicked "best" ones. Second, the noise has been approximated 




(1) 
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by an uncorrelated Bernoulli noise while it is believed that one should either use 
an exponentially correlated dichotomic noise and then take a certain double limit to 
simulate a GWN or directly use good generators of a GWN. We will see that 
the problem of simulating the system described by the equation (JTJ) is unexpectedly 
tricky: Even though this system almost surely never diverges, we will show that its 
discretizations realized by several popular and well-established methods do diverge. 

This paper is organized as follows: In Section |21 we construct stationary probability 
distributions corresponding to the equation ((TJ). Several interacting species may not 
converge to a stable fixed point even if their trajectories are bounded. They may 
display oscillations instead and may not have a stationary distribution. However, if 
there is only one species, finding its stationary distribution is a straightforward task. 
We will show that systems whose deterministic dynamics are completely different may 
possess similar stationary distribution when perturbed by the noise. Then in Sections El 
and E| we will discuss how to simulate the system (JTJ). In particular, we will show that 
there is an "optimal" noise level that maximizes the average lifespan of a discretized 
system. Conclusions are given in Section |3J 



2. The stationary distribution 

The equation (Q) leads to the following Fokker-Planck equation: 
dP(x,t) 



dt 



d a 2 d 2 

-r-[x(l±x)P(x,t)] + T — 



[x 4 P{x,t)] . 



(2) 



A stationary solution to this equation can be easily found by standard methods [SHE]. 
It is, however, instructive to take a different approach. We make a substitution 



1 

y = - 

x 



(3) 



Note that because x is nonnegative almost surely, the substitution (jSJ) is always valid, 
i.e. we do not risk a division by zero almost surely. As a result of this substitution, the 
Fokker-Planck equation (j2J) in the Ito interpretation is transformed into j3] 



dP(y,t) 
dt 



d_ 

dy 



ry ± r 



O' 

v 



P(y,t) 



+ 



Q- 2 d 2 P(y,t) 
2 dy 2 



which, in turn, corresponds to the following Langevin equation: 
/ c 2 ' 

y 



ry ± r + a^(t) . 

y 



(4) 



(5) 



The deterministic part of this equation describes an overdamped motion in a "potential" 

1 



U(y) = -ry ±ry — a \ny . 



(6) 



We can see that the noise interpreted in the sense of Ito introduces a barrier preventing 
y from crossing zero, or preventing the population x from diverging, as predicted in 
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Ref. PQ. Note that this barrier is absent if the equation is interpreted in the sense 
of Stratonovich, and indeed there is no proof that the population x does not explode 
in this case. Finding the stationary solution to the equation (J3J) is now a matter of a 
simple calculation: 



P st (y) =A/Vexp ? y(j/±2) 



or after transforming back to the original variable, 



"stW = — exp 

3v 



_l±2;c 

" r ~ 

ar 



where r = r/a 2 and the normalization constant equals 

4f3/2 

AT 



exp(f)(l + 2f) (l =F Erf (v^)) =F 2v^ 



(7) 



(8) 



(9) 



and Erf(-) is the error function. It is easy to verify that for any positive value of f, and 
for both choices of the sign, the normalization constant M is positive and finite. The 
distribution (jHJ) is the stationary probability distribution corresponding to the equation 
(J2J with the constraint x > 0. It reaches a maximum at x — 

Observe that the "— " sign in (JTJ and the subsequent equations corresponds to 
a noisy logistic equation; the corresponding deterministic system is certainly stable and 
never diverges. The sign "+" corresponds to a system whose deterministic counterpart 
diverges very rapidly. However, the choice of this sign has very little effect on the shape 
and properties of the distribution (jHJ). It is indeed very surprising that two systems whose 
deterministic behaviours differ so dramatically have very similar stationary distributions 
when driven by a GWN. 



3. Numerical simulations 



When one deals with a Langevin equation in the form 

x = f(x)+g(x)£(t) (10) 

and does not know a stationary distribution corresponding to it, the usual way to 
proceed is to use an appropriate discretization in time, solve such system numerically, 
let it "equilibrate" for a certain time and try to reconstruct the (unknown) stationary 
distribution from the numerical trajectories. If £(£) is a GWN, the celebrated Euler- 
Maruyama scheme is most frequently used for this purpose. In this scheme the numerical 
trajectory is generated by 

x n+ i = x n + hf(x n ) + y/hg(x n )r) n , (11) 

where h is the time step and r) n is a discrete Gaussian white noise: (r) n ) = 0, 
(VnVm) = Sum- The Euler-Maruyama scheme is consistent with the equation (jlUj) in 
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the Ito interpretation in the limit h — > + . This scheme is numerically very cheap but 
its strong order of convergence is rather small, 7 = 1/2 

Simulating the noisy logistic equation, corresponding to a "— " in the equation (JTJ), 
does not pose any particular problem. We will, therefore, focus on simulating the system 
whose deterministic counterpart diverges. If we apply the Euler-Maruyama scheme to 
the equation (JTJ), we can see that the population is prevented from exploding by sudden 
jumps down. The amplitude of these jumps is scaled by the square of the current value 
of the population: the higher the system climbs, the lower it can drop. Large jumps 
up are equally possible, but a long series of jumps up is unlikely. If x n is large, even a 
small (and positive) value of 7] n introduces a significant increase in the population size. 
In the next step, however, the probabilities of jumping up and down are equal and even 
a smaller to the absolute value but negative r] n+ i (r] n+ i < 0, \rj n +i\ < Vn) can undo the 
cumulative effect of several jumps up as it is scaled by a much larger factor. We know 
that the ideal, continuous time system neither diverges, nor becomes negative, but the 
discretization scheme (fTT|) applied to the equation (JTJ may not obey these principles. 
If the current magnitude of the population becomes too large or negative, the model 
(JIJ clearly looses any physical meaning. We say that the numerical trajectory becomes 
unphysical either if x n+ \ < or if x n+ \ > X 3> and stop the simulation when either 
of these situations occurs. Observe that the sequence {x n } does not converge to zero 
without becoming negative: The probability 



becomes negligible for < x n -C 1. On the other hand, the probability that x n+ \ 
becomes negative equals 



a consistency check of the Euler-Maruyama scheme, but is nonzero for any h > and 
any x n > 0. This means that after a certain, perhaps very large but finite, number of 
iterations, the Euler-Maruyama scheme will produce a negative, i.e. unphysical, value. 
Thus, the discretized system has a finite lifespan, i.e. the time after which the system, 
when started from a positive initial value, either becomes negative or exceedingly large. 

Observe that if < x n <C 1, probability ()13|) is very small and increases with 
an increasing value of x n . If x n ^> 1, probability (fTH|) reaches a constant value, 
Prob (x n+ \ <0|x n >0) ~ (1 - Eri(Vhr/a))/2, which may be quite large; for example, 
for h = 2 -16 and r/a = 1, Prob (x n+ \<0 | x n 3>l) ~ 0.498. We can see that it is easier 
for a discretized version of system to cross zero if the current value of the process, or 
the current population, is large than when it is small. This conclusion is well confirmed 
by numerical simulations. 





(13) 



This probability goes to zero in the limit h 



+ , which is nothing more than 
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Figure 1. Average lifespan of system JQ discretized by the Euler-Maruyama scheme 
(|llfl as a function of the noise amplitude, a. The lifespan and a are dimensionless. 
The line connecting the points is meant as a guide for the eye only. 



Events in which x n becomes larger than the threshold value X ^> are very rare. 
If x n becomes large, it is more likely for it to decrease in a sudden jump and start its 
way up again, or even become negative when the simulation stops than gradually climb 
up towards the threshold. 

If the noise intensity, a, is very small, the "stopping effect" of the noise is limited: 
the population can build up beyond the threshold or, more likely, to such a value where 
the probability of the next iterate becoming negative is significant despite the small 
value of a. On the other hand, in the limit a — > oo, the probability (fTT?|) goes to 1/2 
and the system becomes unphysical very rapidly. Therefore, we expect the lifespan of 
a system with a very small or very large noise intensity to be short. There should be 
an "optimal" range of the intensity of the noise in which the noise effectively stops the 
system from climbing too high and at the same time the probability ()13j) has a modest 
value. Lifespans of the system with such a noise intensity should be larger than for 
either small or large noise intensities. 

We have performed numerical simulations to verify this. We have generated 
trajectories of the equation (Q) according to Euler-Maruyama scheme with a time step 
h = 2~ 16 ~ 1.5 ■ 10~ 5 . The noise has been generated by Marsaglia algorithm [H] with 
Mersenne Twister as the underlying uniform generator |5] . The simulations were stopped 
when x n+ i became negative or larger than X = 10 20 ; the latter events were very rare. 
For each value of a, the lifespans were averaged over 1024 realizations of the process. 
Results are presented on Figure [T] As we can see, the average lifespan increases rapidly 
as a increases, reaches a plateau, and then gradually decreases. These results show 
clearly that there is a range of "optimal" noise intensities that maximize the lifespan of 
the discretized process. A precise location of the maximum is not particularly important 
as it may depend on realizations of the noise. 
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Figure 2. Numerical (solid line) and theoretical (broken line) histograms of 
distribution |JSJ corresponding to the deterministically divergent case (the "+" sign). 
Clockwise from top-left a = 1, 2, 16, and 128, respectively. The agreement between 
theoretical and numerical histograms is very good. 

The average lifespan displays a strong dependence on the time step used to perform 
the simulations: The average lifespans increase (decrease) with a decreasing (increasing) 
value of the time step, h, and this increase (decrease) can be quite significant. However, 
with very short time steps, the time needed to perform any realistic simulations gets 
prohibitively large. 

Data gathered during the simulations were also used to draw "experimental" 
histograms of the distribution (jHJ). Specifically, from all realizations with h = 2~ 16 and 
lifespans equal or greater than 64, values of x(t) for 32 ^ t ^ 64 were collected, binned 
in intervals of the size 2~ 10 , averaged over the respective number of realizations and 
compared with histograms obtained from the theoretical distribution. Selected results 
are presented on Figure El The agreement between theory and numerical experiment is 
very good for 1/2 < a < 1024, meaning that the Euler-Maruyama scheme reproduces 
properties of the equation (JTJ) well. Results outside this interval are worse because few 
sufficiently long trajectories were available and the statistics was poor; for large a one 
should also use a smaller bin size. 

4. Higher order schemes 

Sometimes there is a need to use a higher order discretization scheme. Several such 
methods are known in literature and in this Section we will discuss performance of two 
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of them. One is the Heun method ^U] which many people, including the present author, 
have applied on several occasions with a great success. However, when this method 
is used to discretize the equation ((H), it leads to negative values of the population 
very rapidly. The reason for this apparent failure is the fact that the Heun method is 
consistent with the equation (fTT)|) in the Stratonovich, not Ito, interpretation. 

The other method, perhaps less popular among physicists, is the Milstein 
scheme 

x n+1 = x n + hf(x n ) + Vhg(x n )r) n + ^g(x n )g'(x n )(r] 2 n -l) , (14) 

where h, r\ n are as in the equation (|TT|) and g'(x n ) = (dg/dx)\ x=Xn . 

The Milstein scheme differs from Euler-Maruyama only in the last term and this 
term vanishes if the noise is purely additive (g = const). Its computational burden 
depends on how difficult it is to evaluate the derivative. In case of a simple polynomial 
the overall cost does not exceed twice that for the Euler-Maruyama method. 

When this method is applied to the equation we find for the probability of the 
next iterate becoming negative 

Prob(x n+ i<0 | x n >0) = - I Erf I — - 1 + Erf I - " ' I I , (15a) 




where 



2 



A = Ah(a x n — rx n — x n ) (15o) 



provided that A > 0; otherwise x n+ \ cannot become negative. This last condition means 
that if 

rVh+ y^ff 2 + hr(r + Aa 2 ) 
x n ^ x min - 2VJ .^ (lb) 

the next iterate must be positive. For a small time step or for a small noise 
intensity, x m i n can be large. This is good news as this prevents the discretized system 
from becoming unphysical. On the other hand, in the limit of a very strong noise, 
Prob(x n+ i<0 | x n >0) — > Erf(l) ~ 0.843 which is much larger than the corresponding 
value for the Euler-Maruyama scheme. The same happens if x n becomes large which is 
likely if the noise intensity is small. There is a trade off between preventing the system 
from becoming unphysical for small values of the population and a rapid divergence 
once the populations builds up sufficiently large. As values of x n are moderate for most 
of the time, we expect that the Milstein scheme would lead to larger lifespans than 
Euler-Maruyama, but the qualitative picture is much the same as in the latter: the 
lifespans for both small and large noise intensities are small and there is a range of noise 
intensities that maximize the average lifespan. These predictions have been confirmed 
numerically. Simulations were run under the same conditions as in case of the Euler- 
Maruyama scheme. Results are presented on Figure El Indeed, the average lifespans first 
increase, then reach a range of elevated values, and eventually decrease. Note that the 
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Figure 3. Same as Fig.^but for the Milstein scheme (|14J) . The lifespan and a are in 
the same units as on Fig. ^ 



maximal lifespans are more than two times larger than those for the Euler-Maruyama 
scheme. The structure seen on Figure 0] (there are apparently two maxima visible) does 
not necessarily correspond to any realistic properties of the iteration (|14j) applied to the 
equation (0). Rather than that, this structure may reflect a very pronounced stochastic 
variability of the system: With the Milstein scheme, an average lifespan of the order of 
2.5 x 10 2 corresponds to a pool of runs with lifespans ranging from very short (~ 10~ 3 ) 
to very large (> 10 4 ) values. 

5. Conclusions 

It is well known that a discretization can dramatically change properties of a dynamical 
system. The deterministic logistic equation is one classic example: its continuous 
version is perfectly regular but after the discretization, the logistic equation may display 
dynamical chaos. Here we have a system — incidentally, closely related to the logistic 
equation — which, when modelled via the continuous time SDE, is mathematically 
proven to remain positive and bounded almost surely, but after a discretization it has 
a finite lifespan after which it becomes negative or diverges. We have argued that, for 
a given time step, there exists a noise level that is "optimal" in that it maximizes the 
average lifespan of the discretized system. Our numerical results strongly support this 
argument. 

We stress that the stabilization effect is observed only in the Ito interpretation. It is 
usually assumed that it is the physical interpretation of the stochastic force that should 
decide which interpretation of the noise, Ito or Stratonovich, should be used pQ. One 
may argue, however, that a robust noise-induced stabilization should be oblivious to the 
noise interpretation. With this respect, Reference |Tj and the present paper provide an 
example of a system in which the consequences of the two approaches differ. 

In this paper we are lucky to numerically simulate a system whose analytical 
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properties are fully known. In practice, it is systems whose analytical properties are not 
known that are examined by numerical simulations. If the simulations produce divergent 
results, one usually concludes that the lifespan of the system under consideration is finite. 
We have shown that such divergences can be merely an artefact of a discretization of 
a SDE. 

Strangely, the above argument can also be reversed: A continuous time SDE is 
a mathematical idealization of a system that is discrete by its very nature. First, 
populations consist of discrete individuals. A fractional x can be interpreted only 
as a fraction of a certain reference population. With x <C 1, this interpretation 
breaks and so may break a model based on a continuous SDE. Second, it is assumed 
that the stochastic force acts continually. This idealization is fully justified when the 
characteristic time scale of the elementary random events is many orders of magnitude 
smaller that the characteristic time scale of a macroscopic process, which is the case 
for example in the classical diffusion. However, when an SDE is used to describe a 
biological or ecological process, the separation of time scales is less perfect. For instance, 
in a biological process the characteristic time of "random" changes in environmental 
conditions can be of the order of seconds, and the characteristic time of the macroscopic 
process, like changes in a population, can be of the order of days. In such a case the 
characteristic time scales are separated by only five orders of magnitude, much as in 
the numerical examples reported above. As a result, a process that is deterministically 
divergent may remain so even after a stochastic perturbation is applied, even though 
solutions to a corresponding continuous time SDE remain bounded almost surely. 

References 

[1] Mao X, Marion G and Renshaw E 2002 Stochastic Process. Appl. 97, 95 

[2] Turchin P 2001 Oikos 94, 17 

[3] Van Den Brocck C 1982 J. Stat. Phys. 31, 467 

[4] Van Kampcn N G 1987 Stochastic Processes in Physics and Chemistry (Amsterdam: North- 
Holland), chapter VIII 
[5] Risken H 1984 The Fokker-Planck Equation (BcrlimSpringcr) 
[6] Gardiner C W 1993 Handbook of Stochastic Methods (Berlin: Springer) 

[7] Klocden P E and Platen E 1992 Numerical Solution of Stochastic Differential Equations 

(Berlin: Springer); Klocden P E 2002 Milan J. Math. 70, 187 
[8] Marsaglia G and Bray T A 1964 SI AM Review 6, 260; Kindcrman A J and Ramage J G 1976 J. 
Amer. Statist. Assoc. 71, 893; Wieczorkowski R and Zielinski R 1997 Komputerowe generatory 
liczb losowych (Warszawa:WNT) (in Polish) 
[9] Matsumoto M and Nishimura T 1998 A CM Trans, on Modeling and Computer Simulation 8, 3 
[10] Mannella R 2002 Int. J. Mod. Phys. C 13, 1177 and references quoted therein. 
[11] Milstcin G N 1974 Theor. Prob. Appl. 19, 557 



